%pylab
import numpy as np
import pylab as plb
import tifffile
#these functions are used for the image overlays.
def overlay_hue_saturation(overlay,bkgnd,
gama = 1.0,
color_map = plb.cm.jet,
alpha = 0.5):
from skimage import data, color, io, img_as_float
gama_correct = lambda x: x**gama
if len(np.shape(bkgnd)) == 2:
img_color = np.dstack((bkgnd, bkgnd, bkgnd))
else:
img_color = bkgnd
if len(np.shape(overlay)) == 2:
overlay_color = color_map(gama_correct(overlay))[:,:,:3]
mask = overlay>0
else:
overlay_color = overlay
mask = np.sum(overlay.astype(uint64),axis = 2)>0
img_hsv = color.rgb2hsv(img_color)
overlay_hsv = color.rgb2hsv(overlay_color)
img_hsv[..., 0] += mask*overlay_hsv[..., 0] * alpha
img_hsv[..., 1] += mask*overlay_hsv[..., 1] * alpha
img_overlay = color.hsv2rgb(img_hsv)
return(img_overlay)
def colorized_image(input_img,
map_point = 0.5,
color_map = plb.cm.jet,
alpha = 0.2,
gama = 1.0,
gain = 1.0):
from skimage import data, color, io, img_as_float
gama_correct = lambda x: x**gama
img = gain*input_img
img = gama_correct(img)
img[img>255] = 255
img = img.astype(uint8)
bg_color = np.dstack((img, img, img))
mask = img>0
map_img = mask*map_point
overlay_color = color_map(map_img)[:,:,:3]
img_hsv = color.rgb2hsv(bg_color)
overlay_hsv = color.rgb2hsv(overlay_color)
#make hue and saturation of img map to color
img_hsv[..., 0] += mask*overlay_hsv[..., 0] * alpha
img_hsv[..., 1] += mask*overlay_hsv[..., 1] * alpha
img_overlay = color.hsv2rgb(img_hsv)
return(img_overlay)
#load images
bkgnd_img = tifffile.TiffFile('./registration/65G06_brightfield_cuticle.tif').asarray()
b1_stack = tifffile.TiffFile('./registration/combined_stack_b1.tiff').asarray()
b2_stack = tifffile.TiffFile('./registration/combined_stack_b2.tiff').asarray()
b3_stack = tifffile.TiffFile('./registration/combined_stack_b3.tiff').asarray()
i1_stack = tifffile.TiffFile('./registration/combined_stack_i1.tiff').asarray()
i2_stack = tifffile.TiffFile('./registration/combined_stack_i2.tiff').asarray()
iii1_stack = tifffile.TiffFile('./registration/combined_stack_iii1.tiff').asarray()
iii3_stack = tifffile.TiffFile('./registration/combined_stack_iii3.tiff').asarray()
iii24_stack = tifffile.TiffFile('./registration/combined_stack_iii24.tiff').asarray()
hg1_stack = tifffile.TiffFile('./registration/combined_stack_hg1.tiff').asarray()
hg2_stack = tifffile.TiffFile('./registration/combined_stack_hg2.tiff').asarray()
hg3_stack = tifffile.TiffFile('./registration/combined_stack_hg3.tiff').asarray()
hg4_stack = tifffile.TiffFile('./registration/combined_stack_hg4.tiff').asarray()
tpd_stack = tifffile.TiffFile('./registration/combined_stack_tpd.tiff').asarray()
tpv_stack = tifffile.TiffFile('./registration/combined_stack_tpv.tiff').asarray()
pr_stack = tifffile.TiffFile('./registration/combined_stack_pr.tiff').asarray()
#make max projections
b1_img = np.max(b1_stack,axis = 0)
b2_img = np.max(b2_stack,axis = 0)
b3_img = np.max(b3_stack,axis = 0)
i1_img = np.max(i1_stack,axis = 0)
i2_img = np.max(i2_stack,axis = 0)
iii1_img = np.max(iii1_stack,axis = 0)
iii3_img = np.max(iii3_stack,axis = 0)
iii24_img = np.max(iii24_stack,axis = 0)
hg1_img = np.max(hg1_stack,axis = 0)
hg2_img = np.max(hg2_stack,axis = 0)
hg3_img = np.max(hg3_stack,axis = 0)
hg4_img = np.max(hg4_stack,axis = 0)
tpd_img = np.max(tpd_stack,axis = 0)
tpv_img = np.max(tpv_stack,axis = 0)
pr_img = np.max(pr_stack,axis = 0)
#colorize each muscle
b_map = plb.cm.Blues_r
i_map = plb.cm.Purples_r
iii_map = plb.cm.Reds_r
hg_map = plb.cm.Oranges_r
tp_map = plb.cm.Greens_r
pr_map = plb.cm.Greens_r
cvls = {'b1':0.01,'b2':0.2,'b3':0.4,
'i1':0.01,'i2':0.2,
'iii1':0.01,'iii3':0.2,'iii24':0.4,
'hg1':0.01,'hg2':0.2,'hg3':0.4,'hg4':0.6,
'tpd':0.01,'tpv':0.1,
'pr':0.01}
b1_colorized = colorized_image(b1_img,color_map = b_map,map_point = cvls['b1'],alpha = 1.0)
b2_colorized = colorized_image(b2_img,color_map = b_map,map_point = cvls['b2'],alpha = 1.0)
b3_colorized = colorized_image(b3_img,color_map = b_map,map_point = cvls['b3'],alpha = 1.0,gama = 0.65,gain = 40.0)
i1_colorized = colorized_image(i1_img,color_map = i_map,map_point = cvls['i1'],alpha = 1.0)
i2_colorized = colorized_image(i2_img,color_map = i_map,map_point = cvls['i2'],alpha = 1.0,gama = 0.7,gain = 10.0)
iii1_colorized = colorized_image(iii1_img,color_map = iii_map,map_point = cvls['iii1'],alpha = 1.0)
iii3_colorized = colorized_image(iii3_img,color_map = iii_map,map_point = cvls['iii3'],alpha = 1.0)
iii24_colorized = colorized_image(iii24_img,color_map = iii_map,map_point = cvls['iii24'],alpha = 1.0)
hg1_colorized = colorized_image(hg1_img,color_map = hg_map,map_point = cvls['hg1'],alpha = 1.0,gain = 2.0)
hg2_colorized = colorized_image(hg2_img,color_map = hg_map,map_point = cvls['hg2'],alpha = 1.0)
hg3_colorized = colorized_image(hg3_img,color_map = hg_map,map_point = cvls['hg3'],alpha = 1.0)
hg4_colorized = colorized_image(hg4_img,color_map = hg_map,map_point = cvls['hg4'],alpha = 1.0)
tpd_colorized = colorized_image(tpd_img,color_map = tp_map,map_point = cvls['tpd'],alpha = 1.0)
tpv_colorized = colorized_image(tpv_img,color_map = tp_map,map_point = cvls['tpv'],alpha = 1.0)
pr_colorized = colorized_image(pr_img,color_map = pr_map,map_point = cvls['pr'],alpha = 1.0)
#combine the images using addition
img_sum = (tpd_colorized+
tpv_colorized+
b1_colorized+
b2_colorized+
b3_colorized+
i1_colorized+
i2_colorized+
iii1_colorized+
iii3_colorized+
iii24_colorized+
hg1_colorized+
hg2_colorized+
hg3_colorized+
hg4_colorized+
pr_colorized)*255
img_sum[img_sum>255] = 255
img_sum = img_sum.astype(uint8)
imshow(overlay_hue_saturation(img_sum[:,:,:3],bkgnd_img,alpha = 1.0))
display(gcf());close()
def mask_stack(top_img,bottom_img,alpha = 0.5,gain = 15):
from skimage import data, color, io, img_as_float
top_hsv = color.rgb2hsv(top_img)
mask = top_hsv[...,2]/np.max(top_hsv[...,2])
mask *= gain
mask[mask>1] = 1
mask[mask<1] = 0
mask = np.dstack((mask,mask,mask))
inverse = 1-mask*alpha
return top_img*mask + bottom_img*inverse
stacked_muscles = reduce(mask_stack,
[tpd_colorized,
tpv_colorized,
pr_colorized,
b2_colorized,
i1_colorized,
hg4_colorized,
iii24_colorized,
i2_colorized,
hg2_colorized,
hg3_colorized,
iii3_colorized,
iii1_colorized,
b1_colorized,
b3_colorized,
hg1_colorized])
gama = 1.5
gama_correct = lambda x: x**gama
cut_gamma = gama_correct(bkgnd_img)+600
cuticle = (np.dstack((cut_gamma,cut_gamma,cut_gamma))*0.03).astype(uint8)
muscles = (stacked_muscles*0.9*255)
sumimg = cuticle+muscles
sumimg[sumimg>255] = 255
sumimg = sumimg.astype(uint8)
imshow(sumimg)
display(gcf());close()
tifffile.imsave('stacked_muscles.tiff',uint8(stacked_muscles*255))
tifffile.imsave('muscle_viz.tiff',sumimg)
mimgs = {'b1':b1_img,
'b2':b2_img,
'b3':b3_img,
'i1':i1_img,
'i2':i2_img,
'iii1':iii1_img,
'iii3':iii3_img,
'iii4':iii24_img,
'hg1':hg1_img,
'hg2':hg2_img,
'hg3':hg3_img,
'hg4':hg4_img,
'tpv':tpv_img,
'tpd':tpd_img,
'pr':pr_img}
model_data = dict()
#construct model
from skimage import measure
from scipy.interpolate import griddata
figure(figsize(5,5))
imshow(sumimg)
for mkey in mimgs.keys():
img = mimgs[mkey]
contours = [measure.find_contours(img,1)[0]]
for n, contour in enumerate(contours):
clen = len(contour[:,0])
x = griddata(np.arange(clen),contour[:,0],np.linspace(clen,20))
x = x[~isnan(x)]
x = hstack((x,x[0]))
y = griddata(np.arange(clen),contour[:,1],np.linspace(clen,20))
y = y[~isnan(y)]
y = hstack((y,y[0]))
model_data[mkey] = vstack((y,x))
plot(y, x,linewidth=0.5,color = 'w')
gca().set_xbound((0,1024))
gca().set_ybound((0,1024))
display(gcf());close()
class Basis(dict):
def __setitem__(self,key,item):
"""overload the __setitem__ method of dict so the transform and inverse
transform matrices are computed when the basis vectors are changed"""
try:
if key in ['a1','a2']:
dict.__setitem__(self,key,item)
A = np.vstack((self['a1'],self['a2'])).T
A_inv = np.linalg.inv(A)
self['A'] = A
self['A_inv'] = A_inv
else:
dict.__setitem__(self,key,item)
except KeyError:
dict.__setitem__(self,key,item)
class GeometricModel(object):
def __init__(self,lines,basis):
self.lines = lines
self.basis = basis
## put lines in barycentric coords
self.barycentric = dict()
for key in self.lines.keys():
coords = dot(self.basis['A_inv'],(self.lines[key]-self.basis['p'][:,newaxis]))
self.barycentric[key] = coords.T
def coords_from_basis(self,basis):
ret = dict()
for key in self.barycentric.keys():
coords = np.dot(basis['A'],(self.barycentric[key]).T)+basis['p'][:,newaxis]
ret[key] = coords
return(ret)
class ModelView(object):
def __init__(self,model):
self.model = model
def plot(self,basis,**kwargs):
lines = self.model.coords_from_basis(basis)
plot_args = {}
plot_args['plot_frame'] = kwargs.pop('plot_frame',True)
plot_args['frame_head_width'] = kwargs.pop('frame_head_width',20)
plot_args['contour_color'] = kwargs.pop('contour_color','w')
kwargs['color'] = plot_args['contour_color']
for line in lines.values():
plot(line[0,:],line[1,:], **kwargs)
if plot_args['plot_frame']:
p = basis['p']
a1 = basis['a1']
a2 = basis['a2']
kwargs['color'] = 'g'
kwargs['head_width'] = plot_args['frame_head_width']
arrow(p[0],p[1],a1[0],a1[1],**kwargs)
kwargs['color'] = 'b'
kwargs['head_width'] = plot_args['frame_head_width']
arrow(p[0],p[1],a2[0],a2[1],**kwargs)
figure(figsize = (10,10))
imfile = tifffile.TiffFile('im_overlay.tiff')
sumimg = imfile.asarray()
###add position of large setae
model_data['e1'] = np.array([[ 170.02788104, 326.71685254],
[ 380.98203222, 919.92627014]])
model_data['e2'] = array([[ 172.83333333, 332.83333333],
[ 551.5 , 164.83333333]])
e1 = model_data['e1']
e2 = model_data['e2']
muscles = dict()
for key in model_data.keys():
if not(key in ['e1','e2']):
muscles[key] = model_data[key]
confocal_basis = Basis()
confocal_basis['a1'] = e2[1]-e2[0]
confocal_basis['a2'] = e1[1]-e2[0]
confocal_basis['p'] = e2[0]
thorax = GeometricModel(muscles,confocal_basis)
thorax_view = ModelView(thorax)
imshow(sumimg)
import copy
thorax_view.plot(thorax.basis)
gca().axis('tight')
display(gcf());close()
# define boolean masks for each muscle
b1_mask = b1_img >0
b2_mask = b2_img >0
b3_mask = b3_img >0
i1_mask = i1_img >0
i2_mask = i2_img >0
iii1_mask = iii1_img >0
iii3_mask = iii3_img >0
iii24_mask = iii24_img >0
hg1_mask = hg1_img >0
hg2_mask = hg2_img >0
hg3_mask = hg3_img >0
hg4_mask = hg4_img >0
tpv_mask = tpv_img >0
tpd_mask = tpd_img >0
pr_mask = pr_img > 0
#Masks for each line
mU56_X_J133 = b1_mask | b2_mask | b3_mask \
| i1_mask | i2_mask \
| iii1_mask | iii24_mask |iii3_mask \
| hg1_mask | hg2_mask | hg3_mask | hg4_mask \
| tpd_mask | tpv_mask \
| pr_mask
mU82_X_J18 = b1_mask | b2_mask | b3_mask \
| i1_mask | i2_mask \
| iii1_mask |iii3_mask \
| hg2_mask | hg3_mask | hg4_mask \
mU82_X_J100 = iii24_mask | hg2_mask
mU82_X_J22 = b1_mask | b2_mask | b3_mask \
|iii3_mask \
| hg1_mask | hg2_mask | hg3_mask | hg4_mask \
mU82_X_J23 = b1_mask | b2_mask | b3_mask \
| iii1_mask |iii3_mask \
| hg2_mask | hg4_mask \
mU82_X_J121 = tpd_mask | tpv_mask
def blur_by_depth(stack):
from scipy.ndimage.filters import gaussian_filter
wf = lambda x: (1.0/(x+30))+1.0
return np.array([gaussian_filter(stack[i],i/2.0 + 3.0)*wf(i) for i in range(shape(stack)[0])])
def model_muscle(stack):
return sum(blur_by_depth((stack>0).astype(float)),axis = 0)
b1_model = model_muscle(b1_stack)
b2_model = model_muscle(b2_stack)
b3_model = model_muscle(b3_stack)
i1_model = model_muscle(i1_stack)
i2_model = model_muscle(i2_stack)
iii1_model = model_muscle(iii1_stack)
iii3_model = model_muscle(iii3_stack)
iii24_model = model_muscle(iii24_stack)
hg1_model = model_muscle(hg1_stack)
hg2_model = model_muscle(hg2_stack)
hg3_model = model_muscle(hg3_stack)
hg4_model = model_muscle(hg4_stack)
tpd_model = model_muscle(tpd_stack)
tpv_model = model_muscle(tpv_stack)
pr_model = model_muscle(pr_stack)
\(X\) is an matrix of \(m\) pixels by \(n\) frames constructed from the signals from \(k\) muscles so that: \(X = WB\) where \(W\) is a \(m\) by \(k\) mixing matrix and \(B\) is a \(k\) by \(n\) matrix of muscle fluorescence over time. We find the best fit to \(B\) as \(W^\dagger X = B\).
import cv2
from scipy.ndimage.morphology import binary_dilation
bk = binary_dilation(mU56_X_J133,iterations = 60).astype(float)
from scipy.ndimage.filters import gaussian_filter
bk = gaussian_filter(bk,50)
mlist = [b1_model,b2_model,b3_model,
i1_model,i2_model,
iii1_model,iii3_model,iii24_model,
hg1_model,hg2_model,hg3_model,hg4_model,
tpd_model,tpv_model,pr_model]#,bk]
model_mtrx = [cv2.pyrDown(cv2.pyrDown(ar)) for ar in mlist]
W = np.array([m.ravel() for m in model_mtrx]).T
W_ = W#= W/np.max(W,axis = 0)
imshow(sum(W_.reshape(256,256,-1),axis = 2),cmap = cm.gray)
display(gcf());close()
#Groups for screen
import flylib as flb
import db_access as dba
import cv2
import cPickle
fly_db = dba.get_db()
U82_X_J121 = flb.Squadron(fly_db,[271,272,273,274,275,276])
U82_X_J18 = flb.Squadron(fly_db,[265,267,268,269,270])##266
U82_X_J100 = flb.Squadron(fly_db,[259,260,261,262,263,264])
U82_X_J22 = flb.Squadron(fly_db,[255,256,257,258])
U82_X_J23 = flb.Squadron(fly_db,[278,279,280,281,282]) ##277
U56_X_J133 = flb.Squadron(fly_db,[244,243,242,241,240,239,238,237])
swarm_dict = {'U56_X_J133':U56_X_J133,
'U82_X_J18':U82_X_J18,
'U82_X_J100':U82_X_J100,
'U82_X_J22':U82_X_J22,
'U82_X_J23':U82_X_J23,
'U82_X_J121':U82_X_J121}
import h5py
testfly = swarm_dict['U56_X_J133'].flies[4]
f = h5py.File(testfly.fly_path + "warped_imgs.hdf5", "r")
X_warped = np.array(f['warped_imgs'])
X = X_warped.T.astype(float)
print shape(X)
print shape(W_)
print numpy.linalg.matrix_rank(W_)
W_pin = numpy.linalg.pinv(W_)
B = dot(W_pin,X)
lns = plot(B.T)
gca().set_xbound(500,2500)
display(gcf());close()
mlist = [b1_model,b2_model,b3_model,
i1_model,i2_model,
iii1_model,iii3_model,iii24_model,
hg1_model,hg2_model,hg3_model,hg4_model,
tpd_model,tpv_model,pr_model,bk]
model_mtrx = [cv2.pyrDown(cv2.pyrDown(ar)) for ar in mlist]
W = np.array([m.ravel() for m in model_mtrx]).T
W_ = W
W_pin = numpy.linalg.pinv(W_)
B = dot(W_pin,X)
lns = plot(B.T)
gca().set_xbound(500,2500)
gca().set_ybound(-0.5,1.5)
display(gcf());close()
X_est = dot(W_,B[:,:5000])
#res = X[:,:5000] - X_est
savetiff = hstack((transpose(X_warped[:5000,:].reshape(5000,256,256),(0,2,1))*5,
transpose(X_est.T.reshape(5000,256,256),(0,2,1))*5))
tifffile.imsave('model_fit4.tiff',transpose(savetiff.astype(uint8),(0,2,1)))
savetiff = hstack((transpose(X_warped[:5000,:].reshape(5000,256,256),(0,2,1))*5,
transpose(X_est.T.reshape(5000,256,256),(0,2,1))*5))
#savetiff += np.min(savetiff)
#savetiff[savetiff > 256] = 256
tifffile.imsave('model_fit5.tiff',transpose(savetiff.astype(uint8),(0,2,1)))
hist(X_est.T[:5000,:].ravel())
max(savetiff.ravel())
for i in range(16):
subplot(16,1,i+1,sharex = gca(),sharey = gca())
plot(B[i,:])
res = X[:,:5000] - X_est
imshow(res[:,0].reshape(256,256),cmap = cm.gray)
import db_access as dba
import flylib as flb
import cv2
import cPickle
output_shape = (1024,1024)
figure(figsize = (8,16))
fly_db = dba.get_db()
swarm = flb.Squadron(fly_db,[244,243,242,241,240,239,238,237])
for i,fly in enumerate(swarm.flies):
#fly = swarm.flies[0]
tiffdata = np.array(fly.fly_record['experiments'].values()[0]['tiff_data']['images'])
pkname = fly.fly_path + '/basis_fits.cpkl'
#print 'here'
f = open(pkname)
img_data = cPickle.load(f)
f.close()
test_basis = Basis()
for key in ['A','p','a1','a2']:
test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
maximg = np.max(tiffdata,axis=0)
col = mod(i,2)
row = i/2
subplot2grid((4,2),(row,col))
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
imshow(cv2.warpAffine(maximg.T,A,output_shape),cmap = cm.gray)
thorax_view.plot(thorax.basis,alpha = 0.5)
gca().set_ybound(0,1024)
gca().set_xbound(0,1024)
gca().set_xticks([])
gca().set_yticks([])
display(gcf());close()
#Groups for screen
U82_X_J121 = flb.Squadron(fly_db,[271,272,273,274,275,276])
U82_X_J18 = flb.Squadron(fly_db,[265,267,268,269,270])##266
U82_X_J100 = flb.Squadron(fly_db,[259,260,261,262,263,264])
U82_X_J22 = flb.Squadron(fly_db,[255,256,257,258])
U82_X_J23 = flb.Squadron(fly_db,[278,279,280,281,282]) ##277
U56_X_J133 = flb.Squadron(fly_db,[244,243,242,241,240,239,238,237])
# define boolean masks for each muscle
b1_mask = b1_img >0
b2_mask = b2_img >0
b3_mask = b3_img >0
i1_mask = i1_img >0
i2_mask = i2_img >0
iii1_mask = iii1_img >0
iii3_mask = iii3_img >0
iii24_mask = iii24_img >0
hg1_mask = hg1_img >0
hg2_mask = hg2_img >0
hg3_mask = hg3_img >0
hg4_mask = hg4_img >0
tpv_mask = tpv_img >0
tpd_mask = tpd_img >0
#Masks for each line
mU56_X_J133 = b1_mask | b2_mask | b3_mask \
| i1_mask | i2_mask \
| iii1_mask | iii24_mask |iii3_mask \
| hg1_mask | hg2_mask | hg3_mask | hg4_mask \
| tpd_mask | tpv_mask
mU82_X_J18 = b1_mask | b2_mask | b3_mask \
| i1_mask | i2_mask \
| iii1_mask |iii3_mask \
| hg2_mask | hg3_mask | hg4_mask \
mU82_X_J100 = iii24_mask | hg2_mask
mU82_X_J22 = b1_mask | b2_mask | b3_mask \
|iii3_mask \
| hg1_mask | hg2_mask | hg3_mask | hg4_mask \
mU82_X_J23 = b1_mask | b2_mask | b3_mask \
| iii1_mask |iii3_mask \
| hg2_mask | hg4_mask \
mU82_X_J121 = tpd_mask | tpv_mask
#for plotting
little_frame = copy.copy(thorax.basis)
little_frame['p'] = thorax.basis['p']/4
little_frame['a1'] = thorax.basis['a1']/4
little_frame['a2'] = thorax.basis['a2']/4
i1_select = i1_mask & ~(b2_mask | b3_mask | b1_mask)
i2_select = i2_mask & ~(iii3_mask | hg4_mask | hg1_mask)
b1_select = b1_mask & ~ (i1_mask | b2_mask | b3_mask)
b2_select = b2_mask & ~ (i1_mask | b1_mask | b3_mask | tpd_mask)
b3_select = b3_mask & ~ (i1_mask | b2_mask | b1_mask | tpd_mask)
iii1_select = iii1_mask & ~ (b2_mask | tpd_mask | tpv_mask | iii24_mask)
iii3_select = iii3_mask & ~ (hg1_mask | hg4_mask | hg3_mask | i2_mask)
iii23_select = iii24_mask & ~ (hg1_mask | hg4_mask | hg3_mask | i2_mask | iii3_mask)
hg1_select = hg1_mask & ~ (hg4_mask | hg3_mask | iii3_mask | iii24_mask)
hg2_select = hg2_mask & ~ (hg4_mask)
hg3_select = hg3_mask & ~ (hg2_mask)
hg4_select = hg4_mask & ~ (iii3_mask | iii24_mask)
tp_select = (tpd_mask | tpv_mask) & ~ (i1_mask | b3_mask | iii1_mask | i2_mask)
figure(figsize(12,4))
for i,mask in enumerate([b1_select,b2_select,b3_select,
i1_select,i2_select,
iii1_select, iii3_select,iii23_select,
hg1_select,hg2_select,hg3_select,hg4_select,
tp_select]):
subplot(2,7,i+1)
imshow(mask)
thorax_view.plot(thorax.basis,lw = 1.0,plot_frame = False,contour_color = 'white',alpha = 0.5)
gca().set_xticklabels([]);gca().set_yticklabels([]);
gca().set_xbound((0,1024));gca().set_ybound((0,1024))
display(gcf());close()
figure()
imshow((hg4_model)+b1_model,cmap=cm.gray)
imshow(sum(model_mtrx,axis = 0),cmap = cm.gray)
plot(B.T)
imshow(X_est.T.reshape(-1,256,256)[0],cmap = cm.gray)
figure()
imshow(sum(res.T.reshape(-1,256,256),axis = 0),cmap = cm.gray)
figure(figsize(8,10))
fly = U56_X_J133.flies[1]
expmnt = fly.experiments.values()[0]
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])
rwing = rad2deg(np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])/5.0)
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])
imgs = tifffile.TiffFile(fly.fly_path + 'warped_imgs.tiff').asarray()
ax = subplot2grid((14,6),(0,0),colspan = 5)
plot(times[:30000],ypos[:30000],rasterized=True)
locator_params(axis = 'y', nbins = 3)
subplot2grid((14,6),(1,0),colspan = 5,sharex = ax)
plot(times[:30000],rwing[:30000],rasterized=True)
locator_params(axis = 'y', nbins = 3)
gca().set_ybound(45,90)
for i,mask in enumerate([b1_select,b2_select,b3_select,
i1_select,i2_select,
iii1_select,iii3_select,iii23_select,
hg1_select,hg2_select,hg3_select,hg4_select]):
subplot2grid((14,6),(2+i,0),colspan = 5,sharex = ax)
mmask = cv2.pyrDown(cv2.pyrDown((mask).astype(float)))>0
plot(times[:30000],mean(imgs[:,mmask],axis = 1),rasterized=True)
if i<11:
gca().set_xticklabels([])
locator_params(axis = 'y', nbins = 3)
subplot2grid((14,6),(2+i,5))
imshow(mask)
thorax_view.plot(thorax.basis,lw = 1.0,plot_frame = False,contour_color = 'white',alpha = 0.5)
gca().set_xticklabels([]);gca().set_yticklabels([]);
gca().set_xbound((0,1024));gca().set_ybound((0,1024))
subplots_adjust(hspace = 0.1)
subplots_adjust(bottom = 0.05)
subplots_adjust(top = 0.99)
subplots_adjust(right = 0.99)
subplots_adjust(left = 0.05)
subplots_adjust(wspace = 0.00)
ax.set_xbound(10,140)
display(gcf());close()
from scipy.optimize import nnls
blist = list()
for i in range(1000):
blist.append(nnls(W,X[:,i]))
B_nn = array([x[0] for x in blist])
X_nn_est = dot(W,B_nn.T)
res_nn = X[:,:1000] - X_nn_est
res = X[:,:1000] - X_est
imshow(X_nn_est.T.reshape(-1,256,256)[0],cmap = cm.gray)
figure()
imshow(sum(res_nn.T.reshape(-1,256,256),axis = 0),cmap = cm.gray)
figure()
imshow(X_est.T.reshape(-1,256,256)[0],cmap = cm.gray)
figure()
imshow(res.T.reshape(-1,256,256)[0],cmap = cm.gray)
plot(B.T)
close()
#this block will load the metadata and warp the images/movies into a common reference frame.
swarm_dict = {'U56_X_J133':U56_X_J133,
'U82_X_J18':U82_X_J18,
'U82_X_J100':U82_X_J100,
'U82_X_J22':U82_X_J22,
'U82_X_J23':U82_X_J23,
'U82_X_J121':U82_X_J121}
def warp_fly_movie(fly):
fname = fly.fly_path + 'meta_data.cpkl'
f = open(fname,'rb');data = cPickle.load(f);f.close()
pkname = fly.fly_path + '/basis_fits.cpkl'
f = open(pkname);img_data = cPickle.load(f);f.close()
test_basis = Basis()
#print data.keys()
for key in ['A','p','a1','a2']:
test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
#compose back to Ca image size
A = vstack((A,[0,0,1]))
S = np.array([[0.25,0,0],[0,0.25,0],[0,0,1]])
Ap = dot(S,A)[:2,:3]
testmovie = np.array(fly.fly_record['experiments'].values()[0]['tiff_data']['images'])
X_warped = [cv2.warpAffine(testmovie[frame].T,Ap,(256,256)).ravel() for frame in range(shape(testmovie)[0])]
X_warped = np.array(X_warped)
import h5py
f = h5py.File(fly.fly_path + "warped_imgs.hdf5", "a")
f['warped_imgs'] = X_warped
f.flush()
f.close()
for swarm in swarm_dict.values():
for fly in swarm.flies:
# os.remove(fly.fly_path + "warped_imgs.hdf5")
print fly.fly_num
warp_fly_movie(fly)
plot(B[:,:5000].T)
X_est = dot(W,B[:,:5000])
print mean(B[:,:5000],axis = 1)
from scipy.optimize import nnls
imshow(Out[51][0])
#this block will load the metadata and warp the images/movies into a common reference frame.
swarm_dict = {'U56_X_J133':U56_X_J133,
'U82_X_J18':U82_X_J18,
'U82_X_J100':U82_X_J100,
'U82_X_J22':U82_X_J22,
'U82_X_J23':U82_X_J23,
'U82_X_J121':U82_X_J121}
img_dict = dict()
output_shape = (256,256)
import cPickle
for swarm_key in swarm_dict.keys():
swarm = swarm_dict[swarm_key]
corr_list = list()
dot_list = list()
max_list = list()
wing_list = list()
norm_list = list()
mean_list = list()
movie_list = list()
rwing_list = list()
ypos_list = list()
for fly in swarm.flies:
try:
#print(fly.fly_num)
fname = fly.fly_path + 'meta_data.cpkl'
#f = open(fname,'rb');data = cPickle.load(f)['cl_stim_metadata'];f.close()
f = open(fname,'rb');data = cPickle.load(f);f.close()
pkname = fly.fly_path + '/basis_fits.cpkl'
f = open(pkname);img_data = cPickle.load(f);f.close()
test_basis = Basis()
#print data.keys()
for key in ['A','p','a1','a2']:
test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
#compose back to Ca image size
A = vstack((A,[0,0,1]))
S = np.array([[0.25,0,0],[0,0.25,0],[0,0,1]])
Ap = dot(S,A)[:2,:3]
#warp the images
dot_list.append(cv2.warpAffine(data['sin_stim_metadata']['dot_img'].T,Ap,output_shape))
corr_list.append(cv2.warpAffine(data['sin_stim_metadata']['corr_img'].T,Ap,output_shape))
max_list.append(cv2.warpAffine(data['sin_stim_metadata']['max_img'].T,Ap,output_shape))
mean_list.append(cv2.warpAffine(data['sin_stim_metadata']['mean_img'].T,Ap,output_shape))
norm_list.append(cv2.warpAffine(data['sin_stim_metadata']['norm_img'].T,Ap,output_shape))
#now warp the movies
new_shape = (shape(data['ave_movie'])[0],output_shape[0],output_shape[1])
warped_ave_movie = zeros(new_shape)
for frame in arange(shape(warped_ave_movie)[0]):
warped_ave_movie[frame,:,:] = cv2.warpAffine(data['ave_movie'][frame,:,:].T,Ap,output_shape)
movie_list.append(warped_ave_movie/mean_list[-1])
rwing_list.append(data['ave_rwing'])
ypos_list.append(data['ave_ypos'])
except Exception as inst:
print inst
img_dict[swarm_key] = {'dot_list':dot_list,
'corr_list':corr_list,
'mean_list':mean_list,
'norm_list':norm_list,
'max_list':max_list,
'movie_list':movie_list,
'rwing_list':rwing_list,
'ypos_list':ypos_list}
ndat = np.array(img_dict[swarm_key]['norm_list'][:]).ravel()
ddat = np.array(img_dict[swarm_key]['dot_list'][:]).ravel()
cdat = np.array(img_dict[swarm_key]['corr_list'][:]).ravel()
norm_max = np.percentile(ndat,95)
dot_max = np.percentile(ddat,95)
corr_max = np.percentile(cdat,95)
norm_min = np.percentile(ndat,0)
dot_min = np.percentile(ddat,5)
corr_min = np.percentile(cdat,5)
def summarize_swarm(swarm_key,swarm_mask,axlist):
axes(axlist[0])
imshow(swarm_mask)
thorax_view.plot(thorax.basis,lw = 1.0,plot_frame = False)
axes(axlist[1])
imshow(np.mean(img_dict[swarm_key]['norm_list'],axis = 0),cmap = cm.gray,vmin = 0, vmax = norm_max)
thorax_view.plot(little_frame,lw = 1.0,plot_frame = False,contour_color = 'orange',alpha = 0.5)
axes(axlist[2])
imshow(np.mean(img_dict[swarm_key]['dot_list'],axis = 0),cmap = cm.gray,vmin = dot_min, vmax = dot_max)
thorax_view.plot(little_frame,lw = 1.0,plot_frame = False,contour_color = 'orange',alpha = 0.5)
axes(axlist[3])
imshow(np.mean(img_dict[swarm_key]['corr_list'],axis = 0),cmap = cm.gray,vmin = corr_min, vmax = corr_max)
thorax_view.plot(little_frame,lw = 1.0,plot_frame = False,contour_color = 'orange',alpha = 0.5)
figure(figsize(8,10))
axlist = [subplot(5,4,i) for i in range(1,5)]
summarize_swarm('U56_X_J133',mU56_X_J133,axlist)
axlist[0].set_ylabel("U56 X J133")
#axlist[1].set_title(r"$\sqrt{diag(XX^\top)}$")
#axlist[2].set_title(r"$Xs^\top$")
#axlist[3].set_title(r"$\frac{Xs^\top}{\sqrt{diag(XX^\top)}||s||}$",y= 1.08)
axlist[1].set_title(r"pixel std")
axlist[2].set_title(r"dot image")
axlist[3].set_title(r"correlation image")
axlist = [subplot(5,4,i) for i in range(5,9)]
summarize_swarm('U82_X_J18',mU82_X_J18,axlist)
axlist[0].set_ylabel("U82 X J18")
axlist = [subplot(5,4,i) for i in range(9,13)]
summarize_swarm('U82_X_J100', mU82_X_J100,axlist)
axlist[0].set_ylabel("U82 X J100")
axlist = [subplot(5,4,i) for i in range(13,17)]
summarize_swarm('U82_X_J22', mU82_X_J22,axlist)
axlist[0].set_ylabel("U82 X J22")
axlist = [subplot(5,4,i) for i in range(17,21)]
summarize_swarm('U82_X_J23', mU82_X_J23,axlist)
axlist[0].set_ylabel("U82 X J23")
for ax in gcf().axes:
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_xticks([]); ax.set_yticks([])
for i in [1,2,3,
5,6,7,
9,10,11,
13,14,15,
17,18,19]:
ax = gcf().axes[i]
ax.set_xbound(15,190);ax.set_ybound(60,236)
for i in [0,4,8,12,16]:
ax = gcf().axes[i]
ax.set_xbound(15*4,190*4);ax.set_ybound(60*4,235*4)
#tight_layout()
display(gcf());#close()
swarm_dict = {'U56_X_J133':U56_X_J133,
'U82_X_J18':U82_X_J18,
'U82_X_J100':U82_X_J100,
'U82_X_J22':U82_X_J22,
'U82_X_J23':U82_X_J23,
'U82_X_J121':U82_X_J121}
m_dict = {'U56_X_J133':mU56_X_J133,
'U82_X_J18':mU82_X_J18,
'U82_X_J100':mU82_X_J100,
'U82_X_J22':mU82_X_J22,
'U82_X_J23':mU82_X_J23,
'U82_X_J121':mU82_X_J121}
def plot_group_sig(sig):
sig_mean = np.mean(sig,axis = 0)
sig_std = np.std(sig,axis = 0)
times = np.linspace(0,2,60*2)
#fill_between(times,sig_mean+sig_std,sig_mean-sig_std,alpha = 0.5)
plot(times,sig.T,color = 'k',alpha = 0.5)#sig_mean)
plot(times,sig_mean,lw = 2.0,alpha = 0.8)
figure(figsize = (8,10))
for j,swarm_key in enumerate(swarm_dict.keys()):
#swarm_key = 'U82_X_J23'
subplot2grid((13,7),(0,j),colspan = 1)
imshow(m_dict[swarm_key])
thorax_view.plot(thorax.basis,lw = 0.5,alpha = 0.5,plot_frame = False)
gca().set_xbound(0,1024);gca().set_ybound(0,1024)
wing_angle = rad2deg(np.array(img_dict[swarm_key]['rwing_list'])/5)
subplot2grid((13,7),(1,j),colspan = 1)
plot_group_sig(wing_angle)
gca().set_ybound(45,75)
def plot_roi(mask):
mmask = cv2.pyrDown(cv2.pyrDown(mask.astype(float)))>0
mdata = np.mean(np.array(img_dict[swarm_key]['movie_list'])[:,:,mmask],axis=2)
plot_group_sig(mdata)
gca().set_ybound(0.9,1.1)
def imshow_roi(mask):
mmask = cv2.pyrDown(cv2.pyrDown(mask.astype(float)))>0
imshow(mmask)
thorax_view.plot(little_frame,lw = 0.5,plot_frame = False,contour_color = 'white',alpha = 0.5)
gca().set_xticklabels([]);gca().set_yticklabels([]);
gca().set_xbound(15,190);gca().set_ybound(60,236)
for i,mask in enumerate([b1_select,
b2_select,
b3_select,
i1_select,
i2_select,
iii1_select,
iii3_select,
hg1_select,
hg2_select,
hg3_select,
hg4_select,
tp_select]):
subplot2grid((14,7),(i+2,j),colspan = 1)
plot_roi(mask)
subplot2grid((14,7),(i+2,6))
imshow_roi(mask)
for ax in gcf().axes:
ax.set_yticklabels([])
ax.set_xticklabels([])
subplots_adjust(hspace = 0.1)
subplots_adjust(bottom = 0.05)
subplots_adjust(top = 0.99)
subplots_adjust(right = 0.99)
subplots_adjust(left = 0.05)
subplots_adjust(wspace = 0.00)
display(gcf())
print(shape(mask_mtrx))
print(numpy.linalg.matrix_rank(mask_mtrx))
imshow(sum(mask_mtrx,axis = 0).reshape(256,256))
X_mean = np.mean(X_warped.astype(float),axis = 0)
#X_Xbar = X_warped/X_mean
X_norm = np.linalg.norm(X_warped.astype(float),axis = 0)
X_Xnorm = (X_warped.astype(float)-X_mean)[~(xwarped == 0)]/X_norm
hist(X_Xnorm[~isnan(X_Xnorm)].ravel())
del(np)
imshow(X_mean.reshape(256,256))
gc.collect()
plot(X_Xbar[500])
imshow(fit_movie.T[9].reshape(256,256),cmap = cm.gray)
tifffile.imsave('fittest2.tiff',fit_movie.T.reshape(10000,256,256).astype(uint8)*10)
tifffile.imsave('fittest_X.tiff',X_warped.reshape(10000,256,256).astype(uint8))
#functions to summarize the sin stimulus protocol
def collect_cycles(ypos,stim_cond,trial_selection_mask,flight_mask):
cycles = list()
trial_selector = zeros_like(ypos).astype(bool)
cycle_selection_mask = ypos>1.5
#trial_selection_mask = around(stim_cond)==2
trial_idxs = flb.idx_by_thresh(trial_selection_mask,0.5)
for trial_idx in trial_idxs:
if not(np.sum((~flight_mask)[trial_idx])):
trial_selector[:] = 0
trial_selector[trial_idx] = 1
cycle_idxs = flb.idx_by_thresh(trial_selector & cycle_selection_mask)
cycles.extend([arange(cycle_idxs[i][0],cycle_idxs[i+1][0])
for i in range(len(cycle_idxs)-1)])
return cycles
def cycle_average(sig,cycle_idxs,times,flight_mask):
"""function extracts the cycle average of the signal given a list of idxs where the cycles
happen, and a mask to exclude non-flight bouts. Resamples so that cycle samples line up.
returns a dictionary of summary stats"""
trial_idxs = [np.arange(x[0],x[0]+210) for x in cycle_idxs]# if (len(x)>0 and (sum(~flight_mask[x[0]:x[0]+200])==0))]
trial_times = [times[idx]-times[idx[0]] for idx in trial_idxs if sum(~flight_mask[idx])==0]
trial_data = [sig[idx] for idx in trial_idxs if sum(~flight_mask[idx])==0]
import scipy.interpolate
fun = scipy.interpolate.griddata
time_resamp = linspace(0,2.,60*2.)
data_resamp = list()
for td,tt in zip(trial_data,trial_times):
f = scipy.interpolate.interp1d(tt, td, kind='linear', axis=0)
data_resamp.append(f(time_resamp))
retdict = {'time_resamp':time_resamp,
'sig_resamp':data_resamp}
return retdict
def wing_correlations(sig,imgs,fly = None):
#calculate some summary images including the correlation image
mean_img = mean(imgs,axis=0)
zeroed_img = imgs - mean_img
max_img = np.max(imgs,axis = 0)
norm_img = numpy.linalg.norm(zeroed_img,axis = 0)
norm_sig = numpy.linalg.norm(sig-np.mean(sig))
dot_img = sum(zeroed_img*(sig[:,newaxis,newaxis]-np.mean(sig)),axis = 0)
corr_img = dot_img/(norm_img*norm_sig)
return {'mean_img':mean_img,
'max_img':max_img,
'norm_img':norm_img,
'dot_img':dot_img,
'corr_img':corr_img}
def calc_fly_stats(fly):
#only did one experiment on these flies
expmnt = fly.experiments.values()[0]
nframes = 30000
#load the movie
movie = np.array(expmnt.exp_record['tiff_data']['images'][:nframes])
flight_mask = np.array(expmnt.exp_record['flight_mask'][:nframes])
#Map some variables
stim_cond = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['DAC2'])[:nframes]
lwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph0'])[:nframes]
rwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])[:nframes]
lmr = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph2'])[:nframes]
aux = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph3'])[:nframes]
xpos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Xpos'])[:nframes]
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])[:nframes]
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])[:nframes]
#idxs associated with sin type stimulus
trial_selection_mask = stim_cond>0
cycle_idxs = collect_cycles(ypos,stim_cond,trial_selection_mask,flight_mask)
#concatinate the stimuli idxs
sin_idxs = hstack(cycle_idxs)
#select the frames from the imaging stack
#that correspond to the sin stimulus
imgs = movie[sin_idxs,:,:]
sig = rwing[sin_idxs]
#correlate with the right wingstroke data
sin_stim_metadata = wing_correlations(sig,imgs,fly)
#now do the same thing for the stripe-fixation epoch at the end of the trial
end_cl_idx = flb.idx_by_thresh(stim_cond<0)[-1]
cl_selection_mask = np.zeros_like(flight_mask)
cl_selection_mask[end_cl_idx] = 1
cl_idxs = squeeze(np.argwhere(cl_selection_mask & flight_mask))
imgs = movie[cl_idxs,:,:]
sig = rwing[cl_idxs]
#again correlate with the right wingstroke data
cl_stim_metadata = wing_correlations(sig,imgs,fly)
#now cycle-average the image stream and left and right wing signals
rdict = cycle_average(movie,cycle_idxs,times,flight_mask)
ave_movie = np.mean(rdict['sig_resamp'],axis = 0)
rdict = cycle_average(rwing,cycle_idxs,times,flight_mask)
ave_rwing = np.mean(rdict['sig_resamp'],axis = 0)
rdict = cycle_average(lwing,cycle_idxs,times,flight_mask)
ave_lwing = np.mean(rdict['sig_resamp'],axis = 0)
rdict = cycle_average(ypos,cycle_idxs,times,flight_mask)
ave_ypos = np.mean(rdict['sig_resamp'],axis = 0)
return {'sin_stim_metadata':sin_stim_metadata,
'cl_stim_metadata':cl_stim_metadata,
'ave_movie':ave_movie,
'ave_rwing':ave_rwing,
'ave_lwing':ave_lwing,
'ave_ypos':ave_ypos}
#calculates and save the 'metadata' for each fly.
#import cPickle
#for swarm in [U82_X_J121, U82_X_J18, U82_X_J100,U82_X_J22, U82_X_J23, U56_X_J133]:
# #for swarm in [U82_X_J18]:
# for fly in swarm.flies:
# try:
# #print(fly.fly_num)
# rdict = calc_fly_stats(fly)
# fname = fly.fly_path + 'meta_data.cpkl'
# f = open(fname,'wb')
# cPickle.dump(rdict,f)
# f.close()
# del(rdict)
# except Exception as inst:
# print inst
for swarm_key in swarm_dict.keys():
print swarm_key
group_movie = np.mean(img_dict[swarm_key]['movie_list'],axis = 0)
movie_std = np.std(group_movie[~isnan(group_movie)])
movie_min = 1-movie_std*8
movie_max = 1+movie_std*8
rescaled = (group_movie-movie_min)*255/(movie_std*16)
rescaled[isnan(rescaled)] = 255/2
rescaled[rescaled>255] = 255
rescaled[rescaled<0] = 0
tifffile.imsave(swarm_key + '.tiff',np.uint8(rescaled))
for i,p in enumerate([0,30,60,90]):
subplot(1,4,i+1)
imshow(group_movie[p],cmap = cm.gray,vmin = movie_min,vmax = movie_max)
mov = np.array(img_dict[swarm_key]['movie_list'])
tifffile.imsave('test2.tiff',uint8(mov[1,:,:,:]*50))
imshow(hg3_mask & ~(hg4_mask | i2_mask | hg1_mask))
shape(img_dict[swarm_key]['rwing_list'])
display(gcf())
plot(rwing[25000])
fm = np.array(expmnt.exp_record['flight_mask'])
argwhere(~fm)[0]
imshow(expmnt.exp_record['tiff_data']['images'][int(argwhere(~fm))])
np.mean(np.array(expmnt.exp_record['tiff_data']['images'][int(argwhere(~fm))]),axis = 0)
F = np.mean(np.array(expmnt.exp_record['tiff_data']['images'])[~fm,:,:],axis = 0)
DF = np.array(expmnt.exp_record['tiff_data']['images'])[fm,:,:]
normed = (DF-F)/F
normed[normed<0] = 0
tifffile.imsave('nstream.tiff', uint8(normed[:500]*50))
plot(rwing[fm])
hist(uint8(normed[:500].ravel()*255))
swarm = U82_X_J18
fly = swarm.flies[0]
fname = fly.fly_path + 'meta_data.cpkl'
f = open(fname,'rb');data = cPickle.load(f);f.close()
data.keys()
pkname = fly.fly_path + '/basis_fits.cpkl'
f = open(pkname);img_data = cPickle.load(f);f.close()
test_basis = Basis()
for key in ['A','p','a1','a2']:
test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
#compose back to Ca image size
A = vstack((A,[0,0,1]))
S = np.array([[0.25,0,0],[0,0.25,0],[0,0,1]])
Ap = dot(S,A)[:2,:3]
imshow(data['ave_movie'][0,:,:])
imshow(warped_ave_movie[50,:,:])
imshow(f(np.linspace(0,2,5000)))
tt
ica_filters, ica_sig = st_ica(mask_signal, ncomp = 7, mu = 0.8)
for i in range(shape(ica_filters)[0]):
tmpimg = zeros((256,256))
tmpimg[mask] = ica_filters[i,:]+abs(np.min(ica_filters[i,:]))
figure()
little_frame = copy.copy(thorax.basis)
little_frame['p'] = thorax.basis['p']/4
little_frame['a1'] = thorax.basis['a1']/4
little_frame['a2'] = thorax.basis['a2']/4
thorax_view.plot(little_frame,lw = 0.3)
imshow(tmpimg)
fly =U_X_J.flies[0]
#load the warped movie
warped_movie = tifffile.TiffFile(fly.fly_path + 'warped_imgs.tiff').asarray()
#only did one experiment on these flies
expmnt = fly.experiments.values()[0]
#Map some variables
stim_cond = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['DAC2'])
lwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph0'])
rwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])
lmr = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph2'])
aux = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph3'])
xpos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Xpos'])
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])
flight_mask = aux < 0.2
trial_selection_mask = stim_cond>0
cycle_idxs = collect_cycles(ypos,stim_cond,trial_selection_mask)
nc = 7
posterior_muscle_mask = i2_mask | iii3_mask \
| iii24_mask | hg1_mask | hg2_mask \
| hg3_mask | hg4_mask
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,posterior_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
posterior_muscle_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
plot_comp_img(A_,posterior_muscle_mask)
def runICA(warped_movie,mask,n_components = 2):
"""get the ICA signals,mixing matrix and unmixing matrix given movie,
mask and expected number of components"""
mask = cv2.pyrDown(cv2.pyrDown(np.float64(mask))) >0
mask_signal = warped_movie[:,mask]
mask_signal = np.float64(mask_signal)
from sklearn.decomposition import FastICA, PCA
ica = FastICA(n_components=n_components,max_iter=1000)
S = mask_signal/mask_signal.std(axis=0)
S_ = ica.fit_transform(mask_signal) # Reconstruct signals
A_ = ica.mixing_ # Get estimated mixing matrix
C_ = ica.components_
return S_,A_,C_
def plot_comp_img(A_,mask):
"""just plot a grid of the component images"""
figure(figsize = (8.5,12))
mask = cv2.pyrDown(cv2.pyrDown(mask.astype(float))) > 0
n_components = shape(A_)[1]
n_rows = int(ceil(nc/3.0))
for i in range(n_components):
col = mod(i,3)
row = i/3
subplot2grid((n_rows,3),(row,col))
tst_img = np.zeros((256,256))
tst_img[mask] = A_[:,i] - np.min(A_)
imshow(tst_img,vmin= 0)#,vmax = 100)
little_frame = copy.copy(thorax.basis)
little_frame['p'] = thorax.basis['p']/4
little_frame['a1'] = thorax.basis['a1']/4
little_frame['a2'] = thorax.basis['a2']/4
thorax_view.plot(little_frame,lw = 0.3)
gca().set_xticks([]);gca().set_yticks([])
gca().set_xbound((0,256));gca().set_ybound((0,256))
def plot_ICA(warped_movie,
mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_):
"""plot a summary of the ICA data from the entire experiment using the components
extracted from a small subset"""
n_components = shape(A_)[1]
n_samples = sum(~flight_mask[:30000])
mask = cv2.pyrDown(cv2.pyrDown(mask.astype(float))) > 0
dat = warped_movie[~flight_mask[:30000],:,:][:,mask]
sigs = np.dot(C_,dat.T)
figure(figsize = (8.5,12))
rwing = rad2deg(rwing/5.0)
#plot the rwing amp cycle average
subplot2grid((n_components+2,5),(0,4))
cd = cycle_average(rwing,cycle_idxs,times)
plot(cd['time_resamp'],cd['sig_mean'])
#plot the rwing example trace
subplot2grid((n_components+2,5),(0,1),colspan = 3)
plot(times[~flight_mask[:30000]],rwing[~flight_mask[:30000]])
ax = gca()
ax.set_ybound((50,90))
raw_sig = np.mean(dat,axis=1)
#plot the cycle average of the raw signal
subplot2grid((n_components+2,5),(1,4))
cd = cycle_average(raw_sig,cycle_idxs,times)
plot(cd['time_resamp'],cd['sig_mean'])
gca().set_yticklabels([])
#plot an example trace of the raw signal
subplot2grid((n_components+2,5),(1,1),colspan = 3,sharex = ax)
plot(times[~flight_mask[:30000]],raw_sig)
gca().set_yticklabels([])
ax = gca()
#itterate through the components
for i in range(n_components):
col = 0
row = i
#first an image of the components
subplot2grid((n_components+2,5),(row+2,col))
tst_img = np.zeros((256,256))
tst_img[mask] = A_[:,i] - np.min(A_)
imshow(tst_img,vmin= 0)#,vmax = 100)
little_frame = copy.copy(thorax.basis)
little_frame['p'] = thorax.basis['p']/4
little_frame['a1'] = thorax.basis['a1']/4
little_frame['a2'] = thorax.basis['a2']/4
thorax_view.plot(little_frame,lw = 0.3)
gca().set_xticks([]); gca().set_yticks([])
gca().set_xbound((0,256)); gca().set_ybound((0,256))
#now the cycle average
subplot2grid((n_components+2,5),(row+2,col+4),colspan = 1)
cd = cycle_average(sigs[i,:],cycle_idxs,times)
plot(cd['time_resamp'],cd['sig_mean'])
gca().set_yticklabels([])
#now an example trace
subplot2grid((n_components+2,5),(row+2,col+1),colspan = 3,sharex = ax)
ax = gca()
plot(times[:n_samples],sigs[i,:])
gca().set_yticklabels([])
ax.set_xbound((40,80))
close('all')
figure();plot(ica_sig[1,:]);
import numpy as np
from numpy import dot,argsort,diag,where,dstack,zeros, sqrt,float,real,linalg
#from numpy import *
from numpy.linalg import eigh, inv, eig, norm
from numpy.random import randn, rand
import time
try:
from scipy.stats import skew
_skew_loaded = True
from scipy.linalg import orth
_orth_loaded = True
except:
_orth_loaded = False
_skew_loaded = False
## TODOs:
## [ ] show how much variance is accounted by pcs
## [X] function to reshape back to movie
## [X] behave nicely if no scipy available (can't sort by skewness then)
## [ ] simple GUI (traits?)
## [ ] masks from ICs or PCs
def pca_svd(X):
ndata, ndim = X.shape
Y = X - X.mean(axis=-1)[:, np.newaxis] # remove mean
U,s,Vh = np.linalg.svd(X, full_matrices=False)
Vh = Vh[:ndata] #only makes sense to return the first num_data
return U,Vh,s
def pca1 (X, verbose=False):
"""
Simple principal component analysis
X as Npix by Nt matrix
X should be normalized and centered beforehand
--
returns:
- EV (Nt by Nesq:esq>0), matrix of PC 'signals'
(eigenvalues of temporal covariance matrix). Signals are in columns
- esq, vector of eigenvalues
"""
tick = time.clock()
Nx, Nt = X.shape
Y = X - X.mean(axis=-1)[:, np.newaxis] # remove mean
C = dot(Y, Y.T)/Nt # temporal covariance matrix
es, EV = eigh(C) # eigenvalues, eigenvectors
## take non-negative eigenvalues
non_neg, = where(es>=0)
neg = where(es<0)
if len(neg)>0:
if verbose:
print "pca1: Warning, C have %d negative eigenvalues" %len(neg)
es = es[non_neg]
EV = EV[:,non_neg]
#S = diag(np.sqrt(es))
#whitenmat = dot(EV, inv(S)) # whitenting matrix
ki = argsort(es)[::-1] # indices for sorted eigenvalues
#U = dot(X, whitenmat) # spatial filters
if verbose:
print "PCA finished in %03.2f sec" %(time.clock() - tick)
return EV[:,ki], es[ki]
# note: returns EVs stored in columns, e.g.
# [[x1, x2, ... xn],
# [y1, y2, ... yn],
# ...]
## Whitening
def whiten_mat1(X):
EV,es = pca1(X)
S = diag(np.sqrt(es))
return dot(EV, inv(S))
def st_ica(X, ncomp = 20, mu = 0.3):
"""Spatiotemporal ICA for sequences of images
Input:
X -- list of 2D arrays or 3D array with first axis = time
mu -- spatial vs temporal, mu = 0 -> spatial; mu = 1 -> temporal
"""
data = X
#data = reshape_from_movie(X) # nframes x npixels
sh = X[0].shape
pc_filters, pc_signals, ev = whiten_svd(data, ncomp)
mux = sptemp_concat(pc_filters, pc_signals, mu)
_, W = fastica(mux, whiten=False)
ica_sig = dot(W, pc_signals)
ica_filters = dot(dot(diag(1.0/np.sqrt(ev)), W), pc_filters)
if _skew_loaded:
skewsorted = argsort(skew(ica_sig, axis=1))[::-1]
else:
skewsorted = range(nIC)
return ica_filters[skewsorted], ica_sig[skewsorted]
def transp(m):
"conjugate transpose"
return m.conjugate().transpose()
pow3nonlin = {'g':lambda X: X**3,
'gprime': lambda X: 3*X**2}
pow3nonlinx = {'g':lambda X,args: X**3,
'gprime': lambda X,args: 3*X**2}
def _sym_decorrelate(X):
"W <- W \cdot (W^T \cdot W)^{-1/2}"
a = dot(X, transp(X))
ev, EV = linalg.eigh(a)
return dot(dot(dot(EV, np.diag(1.0/np.sqrt(ev))),
EV.T),
X)
def _ica_symm(X, nIC=None, guess=None,
nonlinfn = pow3nonlin,
termtol = 5e-7, max_iter = 2e3,
verbose=False):
"Simplistic ICA with FastICA algorithm"
nPC, siglen = map(float, X.shape)
nIC = nIC or nPC
guess = guess or np.random.normal(size=(nPC,nPC))
guess = _sym_decorrelate(guess)
B, Bprev = zeros(guess.shape, np.float64), guess
iters, errx = 0,termtol+1
g, gp = pow3nonlin['g'], pow3nonlin['gprime']
while (iters < max_iter) and (errx > termtol):
bdotx = dot(Bprev, X)
gwtx = g(bdotx)
gp_wtx = gp(bdotx)/siglen
B = dot(gwtx, X.T)/siglen - dot(np.diag(gp_wtx.mean(axis=1)), Bprev)
B = _sym_decorrelate(B)
errx = max(abs(abs(np.diag(dot(B, Bprev.T)))-1))
Bprev = np.copy(B)
iters += 1
if verbose:
if iters < max_iter:
print "Success: ICA Converged in %d tries" %iters
else:
print "Fail: reached maximum number of iterations %d reached"%maxiters
return B.real
def whiten_svd(X, ncomp=None):
n,p = map(float, X.shape)
Y = X - X.mean(axis=-1)[:, np.newaxis]
u, s, vh = linalg.svd(X, full_matrices=False)
K = (u/s).T[:ncomp]
return np.dot(K, X), K, s[:ncomp]
def fastica(X, ncomp=None, whiten = True,
algorithm = 'symmetric',
tol = 1e-3, max_iter = 1e3, guess = None):
n,p = map(float, X.shape)
if whiten:
XW, K, _ = whiten_svd(X, ncomp) # whitened and projected data
else:
XW = X.copy()
XW *= np.sqrt(p)
kwargs = {'termtol':tol, 'nonlinfn':pow3nonlin,
'max_iter':max_iter, 'guess':guess}
algorithms = {'symmetric':_ica_symm, 'deflation':None}
fun = algorithms.get(algorithm, 'symmetric')
W = fun(XW, **kwargs)
if whiten:
S = dot(dot(W,K),X)
else:
S = dot(W,X)
return S, W
def fastica_defl(X, nIC=None, guess=None,
nonlinfn = pow3nonlin,
termtol = 5e-7, maxiters = 2e3):
tick = time.clock()
nPC, siglen = X.shape
nIC = nIC or nPC-1
guess = guess or randn(nPC,nIC)
if _orth_loaded:
guess = orth(guess)
B, Bprev = zeros(guess.shape, np.float64), zeros(guess.shape, np.float64)
iters, minAbsCos = 0,0
errvec = []
icc = 0
while icc < nIC:
w = randn(nPC,1) - 0.5
w -= dot(dot(B, transp(B)), w)
w /= norm(w)
wprev = zeros(w.shape)
for i in xrange(long(maxiters) +1):
w -= dot(dot(B, transp(B)), w)
w /= norm(w)
#wprev = w.copy()
if (norm(w-wprev) < termtol) or (norm(w + wprev) < termtol):
B[:,icc] = transp(w)
icc += 1
break
wprev = w.copy()
return B.real, errvec
### Some general utility functions:
### --------------------------------
def gauss_kern(xsize, ysize=None):
""" Returns a normalized 2D gauss kernel for convolutions """
xsize = int(xsize)
ysize = ysize and int(ysize) or xsize
x, y = mgrid[-xsize:xsize+1, -ysize:ysize+1]
g = np.exp(-(x**2/float(xsize) + y**2/float(ysize)))
return g / g.sum()
def gauss_blur(X,size=1.5):
return signal.convolve2d(X,gauss_kern(size),'same')
def reshape_from_movie(mov):
l,w,h = mov.shape
return mov.reshape(l,w*h)
def reshape_to_movie(X,nrows,ncols):
Np, Nt = X.shape
return X.T.reshape(Nt,nrows,ncols)
def sptemp_concat(filters, signals, mu):
if mu == 0:
out= filters # spatial only
elif mu == 1:
out= signals # temporal only
else:
out = np.concatenate(((1-mu)*filters, mu*signals),
axis = 1)
return out / np.sqrt(1-2*mu+2*mu**2)
def DFoF(X):
"""
Delta F / mean F normalization for each pixel
"""
xmean = X.mean(axis=1).reshape(-1,1)
xmean_z = where(xmean==0)
xmean[xmean_z] = 1.0
out = X/xmean - 1
out[xmean_z,:] = 0
return out
def DFoSD(X, tslice = None):
"""
Delta F / S.D.(F) normalization
"""
if not isinstance(tslice, slice):
tslice = tslice and slice(*tslice) or slice(None)
xsd = X[:,tslice].std(axis=1).reshape(-1,1)
xmean = X[:,tslice].mean(axis=1).reshape(-1,1)
return (X-xmean)/xsd
def demean(X):
"""
Remove mean over time from each pixel
Frames are flattened and are in columns
"""
xmean = X.mean(axis=0).reshape(1,-1)
return X - xmean
def winvhalf(X):
"raise matrix to power -1/2"
e, V = eigh(X)
return dot(V, dot(diag((e+0j)**-0.5),transp(V)))
def msqrt(X):
e, V = eigh(X)
return dot(V, dot(diag(((e+0j)**0.5)), transp(V)))
def mpower(M, p):
"""
Matrix exponentiation, works for Hermitian matrices
"""
e,EV = linalg.eigh(M)
return dot(transp(EV),
dot(diag((e+0j)**p), EV))
def mask4overlay(mask,colorind=0):
"""
Put a binary mask in some color channel
and make regions where the mask is False transparent
"""
sh = mask.shape
z = np.zeros(sh)
stack = dstack((z,z,z,np.ones(sh)*mask))
stack[:,:,colorind] = mask
return stack
def flcompose2(f1,f2):
"Compose two functions from left to right"
def _(*args,**kwargs):
return f2(f1(*args,**kwargs))
return _
def flcompose(*funcs):
"Compose a list of functions from left to right"
return reduce(flcompose2, funcs)
nc = 3
#ica_data = warped_movie[~flight_mask[:10000],:,:]
S_,A_,C_ = runICA(ica_data,b3_mask,n_components = nc)
plot_ICA(warped_movie,
b3_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
plot_comp_img(A_,b3_mask)
print mean(A_,axis = 0)
nc = 12
all_muscle_mask = b1_mask | b2_mask \
| b3_mask | i1_mask \
| i2_mask | iii1_mask \
| iii3_mask | iii24_mask \
| hg1_mask | hg2_mask \
| hg3_mask | hg4_mask
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,all_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
all_muscle_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
plot_comp_img(A_,all_muscle_mask)
nc = 4
anterior_muscle_mask = b1_mask | b2_mask \
| b3_mask | i1_mask \
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,anterior_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
anterior_muscle_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
figure()
plot_comp_img(A_,anterior_muscle_mask)
nc = 2
ica_data = warped_movie[~flight_mask[:5000],:,:]
b3_bool_mask = b3_mask & ~(i1_mask | b1_mask | b2_mask)
S_,A_,C_ = runICA(warped_movie,b3_bool_mask,n_components = 2)
plot_ICA(warped_movie,
b3_bool_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
#display(gcf());close()
fly = U_X_J.flies[0]
#load the warped movie
warped_movie = tifffile.TiffFile(fly.fly_path + 'warped_imgs.tiff').asarray()
#only did one experiment on these flies
expmnt = fly.experiments.values()[0]
#Map some variables
stim_cond = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['DAC2'])
lwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph0'])
rwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])
lmr = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph2'])
aux = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph3'])
xpos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Xpos'])
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])
flight_mask = aux<0.0
trial_selection_mask = stim_cond>0
cycle_idxs = collect_cycles(ypos,stim_cond,trial_selection_mask)
nc = 4
anterior_muscle_mask = b1_mask | b2_mask \
| b3_mask | i1_mask \
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,anterior_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
anterior_muscle_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
figure()
#plot_comp_img(A_,anterior_muscle_mask)
nc = 3
#i1_mask = b1_mask | b2_mask \
# | b3_mask | i1_mask \
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,i1_mask,n_components = nc)
plot_ICA(warped_movie,
i1_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
figure()
#plot_comp_img(A_,anterior_muscle_mask)
nc = 3
#i1_mask = b1_mask | b2_mask \
# | b3_mask | i1_mask \
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,b2_mask,n_components = nc)
plot_ICA(warped_movie,
b2_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
figure()
#plot_comp_img(A_,anterior_muscle_mask)
draw()
plot(proj[1,:])
plot(cycle_average(S_[:,0],cycle_idxs,times)['sig_median'])
s1 = sum(warped_movie[:,mask],axis = 1)
i1_bool_mask = i1_mask & ~(b2_mask)
mask = cv2.pyrDown(cv2.pyrDown(np.float64(i1_bool_mask))) >0
s2 = sum(warped_movie[:,mask],axis = 1)
from sklearn.decomposition import FastICA
ica = FastICA(n_components=1)
S = c_[s1,s2]/c_[s1,s2].std(axis = 0)
S_ = ica.fit_transform(c_[s2,s1]) # Reconstruct signals
A_ = ica.mixing_ # Get estimated mixing matrix
C_ = ica.components_
imshow(b3_bool_mask)
plot(flight_mask)
#the anterior cluster
nc = 4
anterior_muscle_mask = b1_mask | b2_mask \
|b3_mask | i1_mask
S_,A_ = runICA(anterior_muscle_mask,n_components = nc)
plot_ICA(anterior_muscle_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
b2_bool_mask = b2_mask
S_,A_ = runICA(b2_bool_mask,n_components = nc)
plot_ICA(b2_bool_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
i1_bool_mask = i1_mask
S_,A_ = runICA(i1_bool_mask,n_components = nc)
plot_ICA(i1_bool_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
b3_bool_mask = b3_mask
S_,A_ = runICA(b3_bool_mask,n_components = nc)
plot_ICA(b3_bool_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
b1_bool_mask = b1_mask
S_,A_ = runICA(b1_bool_mask,n_components = nc)
plot_ICA(b1_bool_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
iii3_bool_mask = iii3_mask
S_,A_ = runICA(iii3_bool_mask,n_components = nc)
plot_ICA(iii3_bool_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
iii24_bool_mask = iii24_mask
S_,A_ = runICA(iii24_bool_mask,n_components = nc)
plot_ICA(iii24_bool_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
hg1_bool_mask = hg1_mask
S_,A_ = runICA(hg1_bool_mask,n_components = nc)
plot_ICA(hg1_bool_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
hg2_bool_mask = hg2_mask
S_,A_ = runICA(hg2_bool_mask,n_components = nc)
plot_ICA(hg2_bool_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
hg3_bool_mask = hg3_mask
S_,A_ = runICA(hg3_bool_mask,n_components = nc)
plot_ICA(hg3_bool_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
hg4_bool_mask = hg4_mask
S_,A_ = runICA(hg4_bool_mask,n_components = nc)
plot_ICA(hg4_bool_mask,S_,A_,n_components = nc);display(gcf());close()
nc = 3
b2_only_boolmask = i1_mask & ~(b2_mask | b1_mask |b3_mask)
S_,A_ = runICA(b2_only_boolmask,n_components = nc)
plot_ICA(b2_only_boolmask,S_,A_,n_components = nc);display(gcf());close()
### function to itterate through a large list of flies
### and make figures of the fits - probably should migrate this to
### an independent script
import db_access as dba
import flylib as flb
import cv2
import cPickle
fly_db = dba.get_db()
def plot_maxfits(flylist):
numfigs = np.ceil(len(flylist)/6.)
chunks = np.array_split(flylist,numfigs)
for chunk in chunks:
fit_swarm = flb.Squadron(fly_db,chunk)
figure(figsize = (8.5,11))
for i,fly in enumerate(fit_swarm.flies):
output_shape = (1024,1024)
#fly = swarm.flies[0]
tiffdata = np.array(fly.fly_record['experiments'].values()[0]['tiff_data']['images'])
pkname = fly.fly_path + '/basis_fits.cpkl'
f = open(pkname)
img_data = cPickle.load(f)
f.close()
test_basis = Basis()
for key in ['A','p','a1','a2']:
test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
maximg = np.max(tiffdata,axis=0)
col = mod(i,2)
row = i/2
subplot2grid((3,2),(row,col))
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
imshow(cv2.warpAffine(maximg.T,A,output_shape),cmap = cm.gray)
thorax_view.plot(thorax.basis,alpha = 0.5)
gca().set_ybound(0,1024)
gca().set_xbound(0,1024)
gca().set_xticks([])
gca().set_yticks([])
gca().set_title(('Fly%04d')%(fly.fly_num))